function diff = q_1_e(m,psi,param)

    if param.zeta==0
        diff = zeros(size(m));
    else
        f = @(x) ( -(1-param.alpha).*2./param.sigma^2.*param.zeta./x./q_0_e(x,psi,param));
        diff = nan(size(m));
        for i=1:length(m)
            diff(i)= q_0_e(m(i),psi,param).*integral(f,param.m_e,m(i));
        end
    end

end